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ABSTRACT 



Context. Sodium laser guide stars (LGS) are about to enter a new range of laser powers. Previous theoretical and numerical methods 
are inadequate for accurate computations of the return flux and hence for the design of the next-generation LGS systems. 
Aims. We numerically optimize the cw (continuous wave) laser format, in particular the light polarization and spectrum. 
Methods. Using Bloch equations, we simulate the mesospheric sodium atoms, including Doppler broadening, saturation, coUisional 
relaxation, Larmor precession, and recoil, taking into account all 24 sodium hyperfine states and on the order of 100 velocity groups. 
Results. LGS return flux is limited by 'three evils': Larmor precession due to the geomagnetic field, atomic recoil due to radiation 
pressure, and transition saturation. We study their impacts and show that the return flux can be boosted by repumping (simultaneous 
excitation of the sodium Dia and D2b lines with 10-20 % of the laser power in the latter). 

Conclusions. We strongly recommend the use of circularly polarized lasers and repumping. As a rule of thumb, the bandwidth of 
laser radiation in MHz (at each line) should approximately equal the launched laser power in Watts divided by six, assuming a 
diffraction-limited spot size. 

Key words, sodium laser guide stars - resonance fluorescence - Bloch equations - atomic physics 



1. Introduction 

• Laser guide stars (LGS) are becoming essential in providing 
■ artificial beacons for adaptive optics (AO) in large telescopes 
. (Hubin 2009). The generation of 8-10 m class telescopes such 
' as the Very Large Telescope (VLT) on Cerro Paranal, Chile, or 
" the Keck telescopes on Mauna Kea, has been retrofitted with 
LGS and is still operated in many observing programs with nat- 
ural guide stars or no AO at all. The upcoming 30-1- m telescopes 
such as the 42-m European Extremely Large Telescope (E-ELT) 
or the Thirty Meter Telescope (TMT), by contrast, are being de- 
' signed from the start as adaptive telescopes and will require LGS 
in nearly all of their science operation. At the same time, sensing 
subaperture sizes of 20-50 cm and near-kHz AO frame rates, de- 
signed to provide high Strehl ratios in the infrared or even in the 
visible, require unprecedented LGS brightness and hence laser 
power. 

Sodium (Na) LGS at 589 nm are most commonly used be- 
cause of the large fluorescence cross section-abundance product 
of ca. 10"" cm^ x 4x lO'-'m"^, a fluorescence wavelength in 
the visible, and the high altitude of the sodium layer around 80- 
100 km compared e.g. to Rayleigh LGS, (Happer 1994, Table 1), 
which allows one to sense a large fraction of the turbulent at- 
mosphere column above the telescope with a small number 
of guide stars. Unfortunately, powerful diffraction-limited laser 
beams at 589 nm are quite expensive to produce due to the lack 
of solid-state materials that amplify 589 nm or 1178nm light. 
Recently, ESO has demonstrated a frequency-doubled narrow- 
band Raman fiber laser emitting 25 W at 589 nm (Taylor 2009), 



and we expect this technology to significantly improve the ex- 
perimental situation. 

Careful dimensioning of the required laser power and op- 
timization of the laser format (spectrum, polarization, spot size) 
are needed. Numerical simulations are necessary since the exper- 
imental situation so far is unsatisfactory: The sodium layer and 
atmospheric parameters often fluctuate rapidly (Thomas 2008, 
Pfrommer 2009), few reliable powerful lasers at 589 nm have 
been available up to now, and LGS sky experiments lack com- 
monly agreed measurement standards. Sodium cells cannot eas- 
ily simulate mesospheric conditions. This work provides opti- 
mization rules for the case of continuous wave (cw) lasers over 
a wide range of laser powers. 

Sodium LGS take advantage of the 3%' i/2-3^P3/2 dipole tran- 
sition, known as the D2 line. The sodium -S ground state consists 
of two hyperfine multiplets with 8 magnetic substates combined, 
separated by 1.772 GHz, splitting the D2 line into the D2a and 
D2b transition groups, corresponding respectively to the F = I 
and F = 2 Na ground states, where F is the total atomic angular 
momentum quantum number The four multiplets (F -Q . . .3) 
are separated by only 16, 34, and 60 MHz, respectively (16 mag- 
netic substates, see Ungar 1 989). At mesospheric temperatures 
near 190 K, the D2a and D2b lines are Doppler broadened to 
about 1 GHz each, giving rise to the characteristic double-hump 
absorption profile (see for instance Bradley (1992), Fig. 11, or 
Fig. 1 of this work). When the D2a transition is excited by cir- 
cularly polarized light at high irradiance (i.e. optical power den- 
sity, unit W/m^), a large fraction of atoms is pumped into the 
{S ,F -2,m- +2) substate, and the atoms cycle on the transition 
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to the (P, F' -3,m' - ±3) substate, with m the magnetic quantum 
number, so that sodium effectively becomes a two-level system. 
This situation is twofold desirable because of the large transi- 
tion cross section and the directional return light peak towards 
the laser emitter and the telescope. Numerous studies of opti- 
cal pumping have been conducted, e.g. by Happer (1972) and 
McClelland (1*^)8?), and the process is well understood. 

The effective sodium return flux depends on the environmen- 
tal conditions such as collisions with constituent gases, temper- 
ature, and the geomagnetic field. At higher irradiance, the atoms 
are driven away from thermal equilibrium, and it is necessary to 
take into account effects such as saturation, optical pumping, and 
recoil caused by radiation pressure. Understanding the compli- 
cated interplay of these eff'ects and obtaining quantitative values 
of the fluorescence efficiency requires numerical simulations. A 
commonly used method is the solution of the density matrix evo- 
lution (Bloch equation) of a multilevel atom. 

Milonni and Thode ( i / _) simplified the D2 scheme to a 2- 
level Bloch model which they solve in time domain. Bradley 
(1992) simulated the full 24-state density matrix, exciting the 
sodium by a train of short (nanosecond range) laser pulses like 
Milonni, using Runge-Kutta integration for one pulse period 
and exploiting the periodicity. Linear and circular light polar- 
izations were treated. However, both works neglect the geomag- 
netic field. Morris (1994) studied frequency-modulated pulses 
over a wide range of pulse durations and linewidths up to 3 GHz, 
hence spanning the entire Doppler broadened D2 line, employing 
time domain integration, for both linearly and circularly polar- 
ized light. In the end of his paper, he briefly estimates the impact 
of the geomagnetic field, but he does not include Larmor preces- 
sion terms into his Bloch equations. 

To our knowledge, Milonni (1998, 1999) has pubHshed 
the most advanced and detailed Bloch-equation simulation of 
sodium LGS to date, later generalized by Telle (2006, 2008). 
In his 1998 article, Milonni treats the cases of laser pulses that 
are short, comparable, and long compared to the lifetime 
of T = 16.24 ns, using numerical solution methods similar to 
Bradley's. His 1999 publication deals with cw excitation only 
and introduces spin relaxation and Larmor precession into the 
Bloch equations for the first time. He shows that the Larmor 
terms due to the geomagnetic field tend to redistribute the mag- 
netic sublevel populations and hence impede optical pumping. 
All of the works cited in this and the previous paragraph solve 
the Bloch equations separately for a number of different sodium 
velocity classes (100-400) and then perform a weighted average 
over the results, and all neglect recoil. 

An alternative method of simulating atomic fluorescence is 
to use rate equations, either implemented as a set of differen- 
tial equations, (Pique -mU .). (Hiflman 2008), or by employ- 
ing Monte Carlo rate equation techniques, (Kibblewhite 2007a, 
2008). In the latter, one follows a more intuitive approach and 
tracks the time evolution of a single atom, executing photon 
absorption and emission events according to steady-state cross 
sections and branching ratios using a random number genera- 
tor, while accounting for elapsed physical time and recording 
the atomic state history. After a sufficient time span At has been 
simulated, one divides the total number of spontaneously emit- 
ted photons by At to obtain the fluorescent flux. The advantage of 
this method is that it is easy to understand, and complex phys- 
ical events such as collisions with different kinds of particles 
with or without spin exchange can be modeled in a transpar- 
ent way. Also, the transit of optically pumped atoms across the 
Doppler spectrum due to velocity-changing collisions and recoil 
can be easily modeled, hence different velocity classes can be 



properly coupled, including their atomic polarization exchange 
(using Bloch equations, such coupling requires the simultaneous 
solution of sets of 24--dimensional equations per velocity class, 
which we describe in this work for the first time in application 
to sodium LGS, to our knowledge). 

A serious disadvantage of rate equations is that they rely on 
steady-state atomic cross sections, limiting the scope to the sim- 
ulation of events slow compared to the sodium transition life- 
time, (Milonni 1992), and the combination of the excitation- 
emission (Rabi) cycling with Larmor precession is hard to im- 
plement correctly. Furthermore, when compared to atomic Bloch 
equations, rate equations neglect the coherences (off-diagonal 
terms of the density matrix that describe transverse atomic po- 
larization), which is problematic when modeling Larmor preces- 
sion or linearly polarized light (the latter even for B = 0). Finally, 
Monte Carlo rate-equations need considerable CPU time in or- 
der to converge. 

There are also approaches that mix aspects of Bloch and 
rate-equation codes, such as BEACON, which has been adapted 
to model two-step sodium excitation for polychromatic LGS, 
(Bellanger 2004). BEACON neglects atomic collisions and re- 
coil. Guillet de Chatellus (2008) reports that it requires on the 
order of 24 h per run on a 2.6 GHz processor, and that the agree- 
ment with a pure rate-equation model can be good in certain 
cases. 

Throughout this paper, we highlight what may be called the 
'three evils' of sodium LGS, ordered by decreasing importance 

1. Larmor precession, 

2. Recoil (radiation pressure), 

3. Transition saturation (stimulated emission), 

and we suggest ways of mitigating them. Larmor precession 
is powerful enough to thwart optical pumping if the angle be- 
tween the beam and the field lines is large (Drummond 2007), 
(Moussaoui 2009a). The average 50-kHz redshift that one in- 
curs per spontaneous emission due to atomic recoil can lead 
to spectral hole burning (depopulation of the respective atoms 
velocity class) in very bright single-frequency LGS (bandwidth 
< 10 MHz, Hiflman (2008)), but it also offers the opportunity 
of 'snowplowing' the sodium population towards higher veloc- 
ities in the laser beam direction, requiring continuous blueshift- 
ing (chirping) of the laser as first suggested by Bradley ( 1 '^>92), 
at the expense of laser complexity. Kibblewhite (2007b) is work- 
ing on an experimental validation of chirping. Ultimately, at high 
spectral irradiance, stimulated emission becomes relevant, lim- 
iting spontaneous emission. Photons from stimulated emission 
are emitted straight into space, hence becoming useless for LGS. 

There is some ambiguity in the sodium LGS community 
about the term 'saturation', for it is sometimes used to describe 
only transition saturation, and sometimes the depopulation of the 
F = 2 upper ground state towards F = I, known as downpump- 
ing. Once an atom is in the F =1 state, it can only be excited if 
the D2b line is pumped as well, which is known as repumping 
(either by allocating 10-20 % of the laser power to the D2b line, 
or by widening a single laser line to > 2 GHz). Downpumping 
in the absence of repumping becomes more severe with increas- 
ing laser irradiance and reduces the return flux long before the 
onset of stimulated emission, hence true transition saturation. 
Spin-exchange collisions can bidirectionally exchange popula- 
tions between the ground states and within them (transitions be- 
tween the {F = 1, m = +1) and {F = 2, m - ±2) ground sub- 
states are particularly strong). Repumping has already been ex- 
perimentally demonstrated to be able to boost the LGS return 
flux by a factor of 1 .6, (Denman 2006b, Telle (2008 finds a factor 
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2.4 through Bloch simulations), using two separate laser beams 
at different frequencies), and we show in this paper that more 
than a factor of 3 can be achieved. 

To date, the most powerful sodium LGS system is installed 
at the Starfire Optical Range (SOR) at the Kirtland Air Force 
Base near Albuquerque, New Mexico, (Denman 2006a), fed by 
a single-frequency (10-kHz linewidth) cw laser that emits SOW 
with circular or linear polarization. Due to its location, the me- 
dian atmospheric seeing at SOR is significantly worse than at 
most sites of large astronomical telescopes, causing enlarged 
LGS spots in the mesosphere, and consequently limiting the irra- 
diance. The upcoming generation of civilian 20-W-class LGS in 
astronomy, such as for the Adaptive Optics Facility of UT4 of the 
VLT (AOF), is therefore expected to venture into unprecedented 
mesospheric irradiance regimes (much higher laser irradiances 
have of course been applied to gas cells, albeit at higher gas pres- 
sures). Next-generation lasers will emit circularly polarized cw 
radiation, use repumping, have a linewidth of a few MHz, and 
are expected to return on the order of 8 x lO^photons/s/m^ on 
the ground. Quantifying and optimizing their return flux vs. the 
laser parameters is one of the principal purposes of this paper 

So far, we have only focused on photon return flux. What is 
really desired when designing AO systems, though, is a bright 
guide star, hence high luminosity concentrated in a small spot 
size. Compared to the uplink laser irradiance in the sodium 
layer, the above mentioned saturation effects spatially broaden 
the LGS return fluorescence distribution by emphasizing the low 
irradiance regions and dimming the peaks. We can show using 
physical-optics simulations, Holzlohner (2008a), that this effect 
increases the instantaneous spot sizes on a wavefront sensor by 
about 0.1", hence it is not negligible. At large angles between 
the laser beam and the geomagnetic field, downpumping can be 
mitigated by repumping. 

In this work, we present a Bloch-equation method that mod- 
els any alkali atom taking into account spontaneous and stimu- 
lated emission, Larmor precession due to the geomagnetic field, 
arbitrary elliptical light polarization, recoil, on the order of 100 
coupled velocity classes with velocity-changing collisions and 
spin exchange, finite atomic dwell time in the beam (atom re- 
placement), arbitrary laser bandwidth, and repumping. We ne- 
glect nonlinear Zeeman shifts and hyperfine coherences since 
we found them to have a small effect on the result. In contrast to 
the above cited Bloch simulation publications, we directly com- 
pute the steady state solution, which is more efficient than time 
domain solutions (a single run takes about 2 s on a modern PC). 
The program is written in Mathematica and based on the public- 
domain Atomic Density Matrix package (Rochester online). 

In order to validate our results, we have also implemented a 
Monte Carlo rate-equation simulation called Exciter in Matlab 
(Holzlohner 2008b,2009), whose results we compare with the 
Bloch-equation method. 

Section 2 describes the Bloch-equation method. Section 3 
gives details about the simulation parameters. Section 4 presents 
the results, suggesting optimal cw sodium LGS formats, and we 
conclude in Section 5. 

2. Bloch equations 

In order to calculate the observed fluorescence from mesospheric 
sodium atoms, the evolution of the atoms is modeled using the 
optical Bloch equations for the atomic density matrix. The den- 
sity matrix describes the statistical state of an ensemble of atoms 
in the state space of the Na D2 transition. In order to account for 
atoms with different Doppler shifts, the density matrix is also 



considered to be a function of atomic velocity along the laser 
beam propagation direction. (An additional degree of freedom is 
included to account for laser line broadening as discussed be- 
low.) The calculation is semiclassical in the sense that while 
the atoms are treated quantum mechanically, the light fields are 
treated classically (the effect of spontaneous decay must be in- 
cluded phenomenologically since we do not apply field quanti- 
zation). Because the density matrix describes all populations of, 
and coherences between, the 24 Zeeman sublevels making up 
the ground and excited states, the calculation describes, in prin- 
ciple, all saturation and mixing effects for essentially arbitrarily 
large optical and magnetic fields. (In practice, certain coherences 
in the system are neghgible under our experimental conditions 
and can be neglected in order to increase the computational effi- 
ciency.) 

In order to perform numerical calculations, the velocity de- 
pendence of the density matrix is discretized to describe an ap- 
propriate number «, .g. of velocity groups, each with a fixed lon- 
gitudinal velocity. Because coherences between atoms with dif- 
ferent velocities can be neglected, the complete density matrix 
p can be thought of as a collection of n,. g separate but coupled 
density matrices, each of dimension 24 x 24. 

The evolution of the density matrix is given by a generaliza- 
tion of the Schrodinger equation: 

-^p= i-[//,p]+A(p)+/3, (1) 
at in 

where H - H{)+ He + Hb is the total Hamiltonian, with Hq the 
Hamiltonian for the unperturbed energy structure of the atom. 
He - -d ■ E the Hamiltonian for the interaction of the elec- 
tric dipole d of the atom with the electric field E of the light, 
Hb - ■ B the Hamiltonian for the interaction of the mag- 
netic moment yu of the atom with the local magnetic field B, 
ft = h/{2n) where /; = 6.626 x lO"-''* Js is Planck's constant, and 
the square brackets denote the commutator The term A in Eq. (1) 
represents phenomenological terms added to account for relax- 
ation processes not described by the Hamiltonian (Budker 2H( !2). 
In our case these relaxation processes include spontaneous decay 
(omitted from the Hamiltonian due to the semiclassical approx- 
imation), collisional spin relaxation ("S-damping") proportional 
to S^p - S ■ (pS) (Happer 1987), and the exit of atoms from the 
light beam due to motion of the atoms and the beam. In addi- 
tion, there are terms included in A to describe changes in atomic 
velocity due to collisions and light-induced recoil, as well as 
an effective relaxation rate that describes dithering of the laser 
phase in order to simulate a finite bandwidth. These terms are de- 
scribed in more detail below. Each relaxation process described 
by A includes a corresponding repopulation process, so that the 
trace over the density matrix for all velocity groups is conserved, 
corresponding to conservation of the total number of atoms. The 
repopulation process describing the entrance of atoms into the 
beam is independent of p and so is written as a separate term 

Velocity-changing collisions are treated as hard collisions in 
which the velocity of the colliding atom is rethermalized in a 
Maxwellian distribution (no speed memory). The internal state 
of the atom is assumed to be unchanged. 

Light-induced recoil is described phenomenologically by 
causing a fraction Vr/Av^.g. of the excited-state atoms in each ve- 
locity group to be transferred upon decay into the next higher ve- 
locity group. Here v, is the recoil velocity and Av,, ,,. is the width 
of the particular velocity group. This model relies on the fact that 
V, - 2.9461 cm/s (equivalent to a Doppler shift of 50.004kHz) 
is much smaller than the typical value of Av,, ^ . 
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In order to simulate a finite bandwidth laser, a form of phase 
dithering is used (frequency or amplitude dithering can also be 
employed). To avoid resorting to a time-domain calculation, the 
dithering is implemented in the spatial domain: density matri- 
ces are written for two 'regions' with light fields that are n out 
of phase with each other, and relaxation terms are included that 
transfer the atoms between the regions (this doubles the size of 
the system of equations). The model is that of a laser beam with 
very fine 'speckles' of different phases. The result is an effective 
laser spectrum of Lorentzian shape with a width proportional to 
the transfer rate between the regions. This method has been ver- 
ified by comparison to a time-domain model (implemented for a 
nuclear-spinless system) in which the light frequency randomly 
changes with a Lorentzian distribution. Identical results from the 
two methods are obtained for the case in which the rate that the 
light frequency changes is faster than the natural decay rate. 

Equation (1) supplies a linear system of differential equa- 
tions for the density matrix elements, known as the optical Bloch 
equations. Thinking of p as a column vector of «, x 24^ den- 
sity matrix elements, the Bloch equations can be written as 
p - Ap + b, where A and b are a matrix and vector, respec- 
tively, that are independent of p. The vector b corresponds to /3 
and A to the rest of the right-hand side of Eq. (1). 

The laser light field has a frequency component tuned near 
the D2 F - 2 ^ F' transition group (Daa), and may have an 
additional 'repump' component tuned near the F = 1 — > 
transition group (D2b). Thus the matrix A has components that 
oscillate at each of these frequencies. Under the rotating wave 
approximation (Corney 1977), the overall optical frequency is 
removed from A. However, the beat frequency between the two 
light-field components remains. This beat frequency can also be 
removed from the Bloch equations in our case: each frequency 
component interacts strongly with one transition group and very 
weakly with the other, so the weak coupling can be neglected for 
each transition. If, in addition, the small magnetic-field-induced 
mixing between the two hyperfine ground states is neglected, 
the beat frequency can be entirely removed from the evolution 
equations. This makes A time-independent for cw light. To find 
the steady-state density matrix, we can set p = and solve the 
linear system Ap - -b. The vectors p and b have 322 elements 
per velocity class (576 if hyperfine states are not neglected), so 
that the sparse hnear equation system has dimension 32,500- 
65,000 in practice. 

To solve the Bloch equations for a particular set of experi- 
mental parameters, we first choose an appropriate set of veloc- 
ity groups. Since the signal is strongly peaked for atoms whose 
Doppler-shifted resonance frequency is near the light frequency, 
we can obtain more accurate results for a given number of ve- 
locity groups if narrower bins are used for resonant atoms, and 
wider for off-resonant. We have two methods for doing this. 

The first method is to choose two fixed bin sizes, one narrow 
and one wide, and the number of narrow bins to cluster near each 
resonance. The wide bins are then used to take up the rest of the 
Doppler distribution. This method is useful when we don't know 
beforehand what the spectrum of the signal in velocity space is. 

If we have an estimate of the spectrum (obtained using the 
first method), we can refine it using the second method, which 
takes advantage of this knowledge. We create a weighting func- 
tion consisting of three terms: a constant term, which tends 
to make equal-sized bins, a term proportional to the spectrum, 
which makes more bins where the signal is large, and a term pro- 
portional to the magnitude of the second derivative of the signal, 
which makes more bins where the signal changes rapidly as a 
function of velocity. The bin sizes are then found by dividing the 



Table 1. Simulation parameters and their standard nominal val- 
ues 



Variable name Symbol Standard value 



Laser parameters 




Launched laser power in air 


P 


20 w 


Mesospheric laser irradiance 


1 


46 W/m^ 


Central D2 vacuum wavelength 


A 


589.159 nm 


Polarization ellipticity angle 


X 


±nlA (circular) 


Laser FWHM linewidth 


A,/ 





Repumping power fraction 


q 


0.12 


Repumping frequency offset 


^fnb 


1.7178 GHz 


Atomic, atmospheric, and mesospheric parameters 


Outer turbulence scale 


Lo 


25 m 


Geomagnetic field in mesosphere 


B 


0.228 G 


One-way transmission at A at zenith 


T 

^ a 


0.84 


Average mesospheric temperature 




185 K 


Sodium centroid altitude (a.s.l.)* 




92 km 


Sodium column density 


Cnb 


4.0xl0"m-2 


Na beam dwell velocity 




38m/s 


Beam atom exchange rate 


Tex 


l/(6.0ms) 


Na-N2 v.c.c." cross section 


C Na-N2 


0.71x10-" cm^ 


Na-02 v.c.c." cross section 


C Na-02 


0.70x10-'" cm^ 


Weighted v.c.c. ^ rate 


Tvcc 


1/(35 A^s) 


Na-02 spin exchange cross sect, at 


0"Na-02 


0.50xl0-''*cm2 


Weighted spin-exchange rate at r^a 


Js 


1/(490 ytiS) 


Launch telescope (LT) parameters 


Zenith angle 




30° 


LT altitude (a.s.l.)' 




2650 m 


LT aperture 


D 


40 cm 


LT beam radius (1/e') 


w 


0.36O = 14.4 cm 


Launched beam rms wavefront error 


WFE 


100 nm « Alb 


Polar angle of B (laser || z) 


e 


njl 


Azimuth of B (laser || z) 





TXjl 



* a.s.l. = above sea level; ' v.c.c. = velocity-changing collision 



integral of the weighting function evenly into the chosen number 
of bins. 

The linear system is solved using the implementation of 
the iterative BiCGSTAB method (stabilized biconjugate gra- 
dient, van der Vorst 1002) built-in to Mathematica. This is a 
Krylov subspace method in which an initial guess is improved by 
minimizing the residual over a subspace with dimension much 
smaller than that of the full system. The rate of convergence 
of the method is increased by pre-multiplication with a block- 
diagonal preconditioner (approximate inverse of A), obtained by 
setting all terms that connect density matrix elements from dif- 
ferent velocity groups to zero, and then inverting the block for 
each velocity group. 

The fluorescent photon flux per solid angle emitted in a given 
direction can be found from the steady-state solution for p as the 
expectation value of a fluorescence operator (Corney 1977). 

3. Simulation Parameters 

3.7. Determination of parameters 

Many physical constants and the atomic level diagram of Na 
have been summarized by Steck (2009) and will not be re- 
peated here. Table 1 lists further simulation parameters, and we 
walk through it to explain some quantities. We work in MKS 
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units, except that we show magnetic field strength in Gauss 
(1 G = lO"'* T) and atomic cross sections in squared centime- 
ters. 

The launched laser power P equals the laser device output 
beam power, diminished by optical losses in the beam train and 
launch telescope (LT). The value of / = 46 W/m~ is the '50% 
power in the bucket' irradiance Ip/2, as explained in the next sub- 
section. We mostly deal with circularly polarized light (the ellip- 
ticity angle x - arctan(e) = +7r/4 denotes respectively LH/RH 
circular and x - ^ linear polarization; n/l - 2y is the polar 
angle in the Poincare sphere, and s is the ellipticity, hence the 
major-minor axes ratio of the polarization ellipse). The repump- 
ing fraction q is the fraction of the total laser power allocated to 
the repump beam: the Dia beam power is (1 - q)P and the power 
in the D2b beam (tuned Afah above the Dia frequency) is qP. 

The geomagnetic-field strength B has a strong impact on 
the return flux. Its value varies considerably over the world and 
can be computed for different mesospheric altitudes using the 
International Geomagnetic Reference Field (IGRF 2005) model. 
Cerro Paranal in northern Chile, the location of the VLT and the 
reference site for this work (24.6°S, 70.4°W), has B ^ 0.23 G 
at 92 km altitude, Mauna Kea (Hawaii) has B a: 0.35 G, and 
the Starfire Optical Range (Albuquerque, New Mexico) has B x 
0.48 G, about twice the field strength at Paranal. 

The atmospheric transmission r„ at 589 nm was measured at 
Paranal to be 0.89 in photometric nights, (Patat 2i ){)4). We use a 
value that is 5 % lower to account for higher aerosol levels. The 
mesospheric temperature, as well as partial gas densities, can be 
derived using the MSISE-90 (1990) model. 

Our values for the sodium layer centroid altitude //Na and 
column abundance Cwa are obtained from studies taking place 
for over 30 years in Sao Paulo, (Simonich 1979, Moussaoui 
2009b), a site whose latitude differs only by one degree from 
the ESO Paranal Observatory. We believe that the sodium layer 
parameter statistics on the seasonal and daily variations are valid 
for Paranal. 

Atomic collisions have a significant effect on the sodium 
states and hence on the LGS return flux. Since mesospheric 
sodium is rarefied (the total mass of global mesospheric sodium 
is about 600 kg), Na-Na collisions are negligible compared to 
Na-N2 and Na-Oa collisions. Most of these collisions are binary 
(collision of two molecules). One important effect of collisions is 
to change the velocity of the atoms, causing diffusion of optically 
pumped atoms in velocity space, sometimes called v-damping 
(Happer 1987). Since the masses of N2 and O2 molecules are 
comparable to that of Na atoms, one can assume that every colli- 
sion completely randomizes their velocity. The collision rate of 
a gas of particle mass Mi with another gas type of particle mass 
M2 and number density /12 is given by (Wright 2004) 

ri2=«2cr,2^^^- + -j, (2) 

where ks = 1.3807 x lO^^^'J/K is Boltzmann's constant and 
o"i2 = ^(f\ + '"2)^ is the collisional cross section with the ef- 
fective particle radii ri, rT. Measuring these radii for velocity- 
changing collisions is difficult; here we just assume the Van-der- 
Waals radii of rNa = 227 pm, tn, - 250 pm, and ro, - 245 pm 
(van den Berg online), (Fishbane 2005). The effect of other gas 
species is negligible. With these numbers, we have, for example, 
o-Na-N2 = 0-72 X lO"'** cm^ and yNa-N, = 3.98 x 10""'cm-^s-i x 

= l/(62.8yL(s) at the sodium centroid («No = 4.0x10'^ cm"-'). 

The sodium layer has a median FWHM thickness of 11 . 1 km 
and its median centroid lies at //Na - 92 km (Moussaoui 2009b). 



The gas pressure decreases exponentially with altitude and the 
collision rate varies across the layer by about one order of mag- 
nitude, as shown by the MSISE-90 (1990) atmospheric model. 
We compute the mean collision rate of Na with N2 and O2, 
based on a table of hn,, no^, and TNa as functions of altitude, 
weighted by the assumed Gaussian sodium density distribu- 
tion. The result, which we will use throughout this work, is 
Tvcc = yNa-N2,02 - 1/(35 yLfs), which is three times higher than 
Milonni's (1999) assumption of 1/(100 yus). Since the sodium 
abundance and layer thickness and altitude are highly variable 
and its profile often deviates significantly from Gaussian, one 
cannot expect a high accuracy in this parameter 

The other important relaxation mechanism beside v- 
damping is spin-exchange relaxation or S-damping (Happer 
I9K 7, Section 13), in particular for Na-02 collisions. Spin re- 
laxation time constants have been measured between rubidium 
and metastable triplet helium (He*, Dmitriev 2008), as well as 
between rubidium and H2, O2, and N2, (Nagengast 1998) and 
sodium and various gases (Ramsey 1964), (Kartoshkin 1998). A 
major difficulty with such measurements is that the overwhelm- 
ing contribution to S-damping of Na in the mesosphere is due 
to collisional spin exchange with O2, however, in gas cells O2 
oxidizes Na quickly and hence this particular cross section is 
hard to determine experimentally. Theoretical calculations of 
the cross section involve Born-Oppenheimer molecular poten- 
tial curves of doublet/quartet surfaces for Na-02, analogous to 
the singlet/triplet curves for Na-Na and have not yet been carried 
out to our knowledge (Happer 2009). 

Kartoshkin (2009) estimates cr^^ ^^ - 0.5 x lO"'"* cm- at 
185 K, based on spin-exchange cross section measurements of 
Na-He* and 02-He. However, only 1/2 of this cross-section 
is effective in our case (Dmitriev 2008, Eq. 3) Setting 0-12 = 
cr^^ Q /2 in Eq. (2), we obtain ys - l/(680//s) at 92km alti- 
tude, which is close to the value js - 1/(640 /is) that Milonni 
(1999) finds through fitting to experiment. Note, however, that 
Milonni's initial guess was cr^^ ^ = 1.0 x 10"'^ cm^, and he 
does not apply the scaling factor of 1 /2. Performing the same 
sodium density weighted averaging over altitude as above, we 
obtain = 1/(490 //s), which will be used throughout this 
work. We will discuss the sensitivity of the Na return flux to 
variations in js and yvcc in the following section. 

The rms lateral velocity Vy = c/yex describes the sodium atom 
exchange into and out of the beam, where d x 23 cm is the me- 
dian FWHM mesospheric speckle diameter (see the following 
section), and yex is the atom exchange rate. We assume that Vy 
consists of four components that we sum in quadrature since in 
general they have uncorrected directions: 

(a) gas diffusion orthogonal to the beam, 

(b) mesospheric wind orthogonal to the beam, 

(c) beam wander caused by atmospheric turbulence, and 

(d) LGS beam slewing due to star tracking. 

Each of these contributions can be estimated: 

(a) The diffusion coefficient of Na in air (mostly N2) can be cal- 
culated using the Chapman-Enskog formula (Cussler 1 997), 
yielding DNa = 8.56 x 10"* cm^s, close to Milonni's assump- 
tion of Dns = 1.0 X lO^cm^/s. The effective lateral diffu- 
sion velocity across the laser speckle is hence ■\/2/3D/d = 
30m/s, where we apply the factor ^/2/3 since we consider 
only the diffusion orthogonal to the beam. Note that 
scales like V7"Na/«N2- At an altitude of about 104 km the 
mean free path of sodium exceeds d and the atomic motion 
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within the beam is no longer diffusive (the mean lateral bal- 
listic velocity is then 370 m/s). 

(b) We use the horizontal wind model HWM07 which com- 
putes zonal and meridional winds at selected mesospheric 
altitudes and times (Drob 2008). We find a median horizon- 
tal wind speed magnitude at Paranal of 20 m/s. Nocturnal 
fluctuations by a factor of 2 are common. 

(c) Physical-optics simulations show an rms beam wander in the 
mesosphere of about 20 cm = 30°), (Holzlohner 2008a), 
a value supported by analytical approximations (Andrews 
2005, Ch. 12, Eq. (51)). The beam wanders on the time scale 
of Tbeam - c//(atmospheric wind speed) x 23 cm / lOm/s - 
23 ms, leading to a beam wander velocity of 8.7 m/s. 

(d) The tracking speed of the laser beam in the mesosphere at 
^ = 30° is about (27r/24h) x L = 7.6 m/s, where L = 
XiHt^a - //tele) ~ ^ 90 km is the line-of-sight distance to 
the guide star centroid with the geometrical length exten- 
sion factor (airmass) X - sec(^), equaling the secant of the 
zenith angle ^. 

Summing these four velocities in quadrature yields Vy - 38 m/s 
and hence -/ex - 1/(6. 0ms). 

3.2. Mesospheric spot size 

Given a certain launched laser power, the spot size and shape 
determine the mesospheric irradiance and thus have a strong im- 
pact on the return flux, particularly in the presence of strong op- 
tical pumping and for small magnetic field polar angles 6. The 
instantaneous mesospheric laser spot size has been simulated us- 
ing physical optics (Holzlohner 2008a). Atmospheric turbulence 
produces a fast changing (timescale a few milliseconds) speckle 
pattern on the sky due to diffraction, where the speckles have a 
FWHM diameter of about 



Table 2. Mesospheric spot parameters for a 40 cm launch tele- 
scope 



2w' 



(3) 



and w is the 1 /e^ irradiance radius of the laser beam at the pro- 
jector. The number of speckles and their beam wander are gov- 
erned by the seeing, and the pattern changes on the time scale of 
a few milliseconds. Note that the medium/long-term spot size as 
observed on the ground is significantly larger than d, and it does, 
in contrast to d, depend directly on the Fried length ro (Fried 
1 966), which is a size scale over which atmospheric phase vari- 
ations remain below 27r. For a typical ro - 5-25 cm, assuming 
an outer turbulence scale of Lq - 25 m, there is a single strong 
central speckle and a couple of much weaker satellite speckles. 

The LGS return flux depends in general nonlinearly on the 
mesospheric laser irradiance /, and the return light has to pass 
the turbulent atmosphere again on the downlink to reach the tele- 
scope; thus the irradiance distribution I{x,y), as it would appear 
on a screen at the mesosphere, cannot be directly observed. To 
account for the effect of turbulence on the return light and char- 
acterize speckle pattern statistically, we collect statistics over 
many simulated realizations of Iix,y) for different turbulence 
phase screens. We rewrite the usual optical power integral 



dl. 



(4) 



as an integral over irradiance up to the maximum value I^^^, 
where A is area orthogonal to the beam and (Tg)^ is the atmo- 
spheric transmission along the line of sight (zenith angle 0. The 



seeing (") 


WFE (nm) 


//>/2(W/m-) FWHMeff(cm) 


d (cm) 


0.6 


70 


79.1 30.0 


22.5 


0.6 


100 


59.1 34.7 


24.3 


1.0 


70 


43.5 40.4 


23.8 


1.0 


100 


33.5 46.1 


26.1 



quantity dA/dl can be thought of as a histogram of mesospheric 
laser irradiance. We define Ip/2 by 



rhii dA r'»"« dA 



(5) 



representing an irradiance of '50 % power in the bucket', which 
will be convenient to compute laser efficiencies later (note that 
usually /max s> 2Ip/2). In order to model a realistic LGS beam 
scenario, we assume a /) = 40cm LT, a launched Gaussian beam 
of l/e- (in irradiance) radius w - 0.72D/2 = 14.4 cm, two 
different seeing values of 0.6" and 1.0" (site monitor seeing at 
500 nm and zenith), and total wavefront errors (WFE) of the 
launched beam of 70 nm rms /1/8) and 100 nm rms (s; A/6). 
We simulate the actual mesospheric irradiance distribution and 
Ip/2 using physical optics, and record the equivalent beam di- 
ameter FWHMeff that a single Gaussian spot of the same meso- 
spheric power and same Ip/2 would have, as well as the speckle 
size d, as in Holzlohner (2008a). Table 2 summarizes the results. 
The mean of Ip/2 on lines 2 and 4 equals 46W/m^, which we 
will use as the reference irradiance in this work. The median 
values of d in the last column are somewhat above the analytical 
prediction of Eq. (3) of 2 1 .2 cm, mostly due to the finite WFE. 

3.3. Flux quantities 

In the following, we define five flux quantities used in this arti- 
cle: 

- Our simulations compute the return flux per solid angle 
as the number of photons per atom and per unit time spon- 
taneously emitted in the direction of the launch telescope, as 
observed in the mesosphere, with units of photons/s/sr/atom. 

- In order to compare different laser formats at similar irradi- 
ances, one can divide ^' by the laser irradiance (in W/m^) 
in the mesosphere to derive the specific atomic return flux 
tf/ = ^'/I in units of photons/s/sr/atom/(W/m-), where area is 
measured in the mesosphere orthogonal to the beam. Most 
of our results are expressed in this quantity. 

- To compute the expected return flux on a receiver at a dis- 
tance L from the fluorescing atoms, we begin by integrating 
'I' over the receiver area A orthogonal to the beam, subtend- 
ing the solid angle A/L- to obtain the quantity F,„. 

- The photon flux on the detector <t) (unit photons/s/m^) at the 
sodium column abundance Cns (unit atoms/m-) equals 



<5 = 



(6) 



The airmass X appears in the numerator because the sodium 
column along the laser beam scales like the airmass (note 
that the distance L grows by the same factor). We do not ac- 
count for laser power depletion with increasing propagation 
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Fig. 1. Absorption cross section cr vs. detuning from the Na 
line center. Solid blue: simulation at / = 0.001 W/m^ {q - 0), 
dashed red: analytical cross section, green: simulation of effec- 
tive cross section at standard conditions (7 = 46 W/m^, ^ = 0.12). 
Inset: reference coordinate system. 



distance in the sodium layer because only about cr x Cns ~ 
10"'^m^x4xl0'-' m"- ^ 4% of the laser photons inter- 
act with any sodium atoms under the standard conditions of 
Table 1, the other 96 % travel unused into space. 
In order to compare LGS systems at similar powers but 
different laser formats, beam spot sizes in the mesosphere, 
atmospheric transmission, and zenith angle, we can di- 
vide F,„ by the laser power in the mesosphere P(Ta)^ 
to arrive at the figure of merit quantity s^e in units of 
photons/s/W/(atoms/m-) 



<DL- 



(7) 



where (Tn)'^ appears quadratically in the denominator of the 
last term, accounting for both up- and downlink. The quan- 
tity See is hence the photon flux on the ground, corrected 
for its dependence on sodium centroid height and abun- 
dance, airmass, atmospheric transmission, and, to first order, 
launched laser power (since Sce depends strongly on theta, 
we will only compute it for 9 - njl in this work). We will 
compare LGS laser beam efficiencies based on Sce- Note that 
in our definition Sce is not a slope efficiency (which would be 
proportional to dF„,ldP); see also d'Orgeville (2000). 

4. Results 

We compute i/r using the method and parameters described 
above. Using the standard conditions of Table 1, we obtain 
(A - 258ph/s/sr/atom/(W/m-). In the following, we first provide 
some numerical validations including scans of (/'(/), followed 
by optimization studies of 4> when varying light polarization 
ellipticity angle repumping frequency offset A/o^, repump- 
ing power fraction q, and laser linewidth A/. Unless otherwise 
noted, the standard conditions of Table 1 apply. 

4.1. Some numerical validations 

Figure 1 shows a comparison between the analytical Na absorp- 
tion cross sections given by Eqs. (10,11) in (Milonni 1998) at 
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Fig. 2. Atomic velocity distribution under standard conditions 
(but 9 - q = Q), blue: Exciter, green: Doppler profile, ma- 
genta: Bloch equations, cyan: return light spectrum in the atomic 
frame. Inset: zoom on line center. 



T - 185 K (dashed red line, the sum of two Gaussians) and the 
simulated effective cross section 



hvW 



(8) 



where v - c/Ais the center frequency of the Daa line, with c the 
vacuum speed of light, and W is the actual rate of spontaneous 
emissions at irradiance I (solid blue curve; the inverse of W is 
known as the cycle time). The abscissa shows detuning of the 
(main) laser line from the D2a line center The agreement is ex- 
cellent. We add the caveat that this agreement is necessary, but 
not sufficient to prove the validity of our algorithm since effects 
like Larmor precession, recoil, and stimulated emission have no 
effect at low irradiance as the Na atom is in thermal equilibrium. 

The steady-state spontaneous emission rate W cannot exceed 
1/(2t) = 1 /(32.5 ns), and hence cr tends to zero in the limit of 
infinite /. At 7 = 46 W/m- and q = 0.12 (green curve), the reduc- 
tion in cr is still modest, however We render the tail of the green 
curve towards the D2b line center at 1.772 GHz dashed, where 
our computational simplification of letting the 'D2a laser line' 
only excite the D2a Na transitions, and analogously allowing the 
'repumping' line to only excite the D2b transitions, whenever 
repumping is used (q > 0), breaks down. Once the detuning ap- 
proaches the D2b line center, this assumption obviously becomes 
invalid. By contrast, the blue curve was computed for ^ = and 
without using this simplification, and it is valid for any frequency 
ofi'set. 

The inset in Fig. 1 sketches the definition of the spherical 
angles and of the B-vector in a coordinate system where the 
laser beam is projected along the z-axis. The major axis of the 
polarization ellipse for non-circular polarization is parallel to x. 

Figure 2 shows the simulated atomic velocity distribution 
under standard conditions, except that Larmor precession and 
repumping are absent (0 = ^ = 7ex = 0). The abscissa shows rel- 
ative atomic velocity away from the receiver in frequency units 
(proportionality constant I /A). The blue line represents the oc- 
cupation histogram from the Monte Carlo rate-equation simula- 
tion Exciter after sampling 10* different atomic velocities within 
the Doppler distribution (8.5 x 10* Monte Carlo steps simu- 
lating 26.6 s of physical time). The green curve is a Gaussian 
with the theoretical Doppler FWHM width of 1.033 GHz at 
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T^Na - 185 K, representing the thermal equilibrium, and the ma- 
genta curve shows the occupation probability obtained from the 
Bloch-equation simulation. Finally, the cyan curve depicts for 
comparison the simulated return flux spectrum in the atomic 
frame, which in the present case of single-frequency excitation 
is close to the sodium natural line shape of a Lorentzian with a 
FWHM of 1 /(Ittt) -9.8 MHz (plotted at arbitrary vertical scale). 

The effect of spectral hole burning is quite striking; in fact 
the occupation at the D2a line center is depleted to 64 % below 
the green Doppler curve. The atomic population, as seen from 
the telescope, is blue-shifted within about one velocity class. We 
observe that spectral hole burning due to recoil is in general ex- 
acerbated if repumping is applied, presumably due to the larger 
number of excitations per time. On the other hand, hole burning 
is mitigated when the laser bandwidth is extended at constant 
laser power because of the reduced spectral irradiance. 

In Exciter, we increment the Doppler shift of the simu- 
lated atom by 50 kHz each time a spontaneous emission occurs, 
which is correct on average. The agreement between Exciter and 
the Bloch equations is excellent, giving us confidence that the 
Gaussian velocity distribution and recoil are properly modeled. 

At / = 46 W/m-, Exciter simulates a spontaneous emission 
every 9.5 yus on average. Of all emissions, 10.2 % are stimulated 
and 89.8 % are spontaneous, and 70.6 % of all atomic transitions 
occur along {F -2,m-2) (F' = 3, m' = 3), indicating efficient 
optical pumping. We conclude that besides suffering from recoil, 
we are also starting to saturate this transition. At / = 100 W/m^, 
17.9 %, and at / = 1000 W/m^, 55.2 % of afl emissions are stim- 
ulated, respectively. On average, one spontaneous emission oc- 
curs every 5.1 yus and 1.1 /is for these respective irradiances. 

Figure 3 shows three semilogarithmic plots of for the 
standard conditions of Table 1 unless noted otherwise, but ne- 
glecting exchange with thermal Na atoms outside the beam 
(Tex = 0, which is not implemented in Exciter, but causes only 
a small difference). We hence simulate a single-frequency laser 
(A/ — > 0) tuned to the peak of the D2a line. Solid curves show 
if/ obtained from the Bloch method, and the symbols indicate the 
results of Exciter. The green and red curves and symbols rep- 
resent the case of laser beam parallel to the geomagnetic field 
{0 - 0) with and without repumping, hence in the absence of 
Larmor precession, respectively. The blue and magenta curves 
and symbols are the same as green and red, respectively, but with 
the laser orthogonal to the field {0 - nil). Figure 3a shows the 
case of circular, and Fig. 3b that of linear polarization, the blue 
and magenta curves for -(f) -njl (note the difference in vertical 
scale). 

We can make a number of interesting observations. First, the 
impact of the magnetic field is profound and reduces the return 
flux strongly, at some irradiances by several times. At very low 
irradiance (/ = 10 ~ W/m^), the atom is in thermal equilibrium 
and all magnetic sublevels are nearly equally populated. Optical 
pumping sets in with increasing / if using circular polarization, 
but Larmor precession is powerful enough to completely sup- 
press it at = ;r/2, as evident from the monotonically falHng 
magenta curve. Conversely, in the absence of Larmor precession 
ij/ strongly peaks near 1-2 W/m-. 

The green and blue curves portray the 'healing' effect of 
repumping {q - 0.12, i.e., 12% of the laser power shifted to 
Dab). In the absence of Larmor precession, the peak is more 
pronounced and shifted to higher / (green). The highest im- 
pact occurs when the influence of the magnetic field is highest 
(9 - njl), where ij/ remains shallow until 80 W/m- and then de- 
cays gently. This finding refutes the often heard notion in the 
LGS community that 'one should always stay away from the cir- 
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Fig. 3. Specific return flux i^(/). Lines: Bloch equations, sym- 
bols: Monte Carlo rate-equation simulation Exciter, a) Green, 
crosses: - Q, q - 0.12; red, diamonds: - Q, q - Q\ blue, 
squares: - njl, ^ = 0.12 (standard conditions); magenta, cir- 
cles: - njl, q - Q (all with circular polarization and yex =0). 
b) Same as a), but for linear polarization, c) Red and magenta 
lines, diamonds and circles: as in a); black, crosses and cyan, 
squares: same, respectively, without recoil. 



cular polarization saturation irradiation (/sat.ciicuiai = 62.6 W/m^, 
Steck (2009)), to avoid saturation' . 

At very low /, repumping is ineffective and even slightly de- 
creases the return flux due to the smaller cross section of the 
Dab transition. At high /, repumping can readily compensate 
the degrading effect of Larmor precession; in its absence, i/r de- 
cays to zero. For linear polarization, this decay even occurs for 
0-Q (red curve) due to depopulation of the upper ground state. 
Conversely, circular polarization 'rescues' many atoms into the 
{S ,F - 2,m - +2) <-» {P, F' - 3, m' - +3) transition, in which 
they are safe from downpumping. 
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Fig. 4. Magnetic field impact, a) Blue curve: i// as a function of 
the Larmor frequency 1/tl at / = 1 W/m^. Red dashed curve (al- 
most obscuring the blue curve): fit (sum of two Lorentzians). 
Blue diamonds: Exciter. Solid green (dashed magenta) lines: 
1/tl at Paranal (SOR). h) F-2 fractional sublevel populations, 
solid (dashed) curves: Bloch (Exciter) simulation. 

Note that besides increasing the return flux by several times, 
a flatter function i^(/) also leads to smaller observed spot sizes 
since spatial power-broadening is reduced, a very welcome 
bonus. 

Another interesting observation is that iff has a local maxi- 
mum near / = 80 W/m^ for the case of = 7t/2 (bump in blue 
and magenta curves). Here we can observe the effect of a mag- 
netic resonance in which the Rabi cycle time in the most strongly 
excited velocity classes becomes similar to the Larmor preces- 
sion time of 



where j-ig is the Bohr magneton and gf is the Lande factor that 
depends on F (gf = +1/2 for the F =1, 2 sodium ground states, 
respectively). On Paranal, ti » 6.3 jUS (159 kHz) and at SOR we 
find tl ^ 3.0 //s (333kHz, B = 0.48 G). On resonance, the pre- 
cessing atoms tend to have the same orientation towards the laser 
whenever the excitations occur (same m), limiting the population 
mixing and consequently mitigating downpumping. We caution 
the reader, however, that the location and height of the bump 
depend on the simulation parameters (e.g. the bump almost van- 
ishes foryvcc = l/(100//s)orB = 0.48 G). 

The agreement between the Bloch calculation and Exciter 
is very good in almost all cases, given the completely different 
nature of the two methods, with the Bloch code usually yield- 
ing the higher values of tfr. Exciter finds the magnetic resonance 
bump at smaller / due to difficulties with the proper simula- 
tion of two concurrent effects on similar, but not equal, time 
scales in our Monte Carlo rate-equation. Also the agreement be- 
tween the green curve and the green crosses is somewhat poor 



Table 3. Sensitivity of to a 1 % perturbation in some simula- 
tion parameters under standard conditions. 



Name of perturbed variable Symbol tff change ( %) 



Mesospheric laser irradiance 


/ 


0.025 


Geomagnetic field in mesosphere 


B 


-0.33 


Average mesospheric temperature 




-0.48 


Beam atom exchange rate 


Tex 


0.0025 


Weighted v.c.c* rate 


Tvcc 


0.15 


Weighted spin-exchange rate at 


7s 


0.013 


Recoil frequency (50 kHz) 




-0.24 



* v.c.c. = velocity-changing collision 



around / = 0.2-20 W/m^, presumably due to the simphfied way 
in which Exciter models S-damping and/or the absence of co- 
herences. 

Figure 3c shows the impact of recoil. The red and magenta 
curves and symbols are as in subplot a), while black and cyan 
crosses and squares denote the same, respectively, but neglecting 
recoil. The magnetic resonance is better visible in the cyan than 
in the magenta curve. Recoil leads to a significant reduction in iff 
above 2 W/m^. 

Figure 4a shows as a function of the Larmor frequency 
1/ti (blue curve) at / = 1 W/nr. The red dashed curve, overlap- 
ping the solid blue curve, is a fit function composed of the sum of 
two Lorentzians of different widths, centered at the origin, plus 
a constant term. The (half-)width of the narrower Lorentzian 
(517 Hz = 0.257s) is determined by the S-damping resonance, 
and that of the broader Lorentzian (17.9 kHz = 0.63')'vcc) by 
the velocity-changing collision rate, and the widths change pro- 
portionally when varying ys or yvcc- However, both resonances 
are somewhat power-broadened. The geomagnetic field is strong 
enough at / = 1 W/m^ to place us on the lowest terrace of the 
blue curve, as indicated by the vertical lines. The blue diamonds 
show the result of Exciter for comparison, also exhibiting the 
terraces. 

Figure 4b shows the corresponding relative populations of 
the five S, F - 2 upper ground state sublevels (m = -2 . . .2), 
where the solid lines indicate Bloch equations and dotted lines 
Exciter For 1/tl > 300 Hz, the sublevel populations collapse 
to the same value due to Larmor-induced sublevel mixing. 
Conversely, for 1/tl< 300 Hz, the populations diverge and the 
laser pumps the m - 2 sublevel (blue) most strongly. When in- 
creasing the irradiance to / = 46 W/nr (not shown in the plot), 
the curves in b) would look similar, but the divergence point 
shifts upward to a; 1 kHz. Furthermore, when plotting the sub- 
level populations at / = 46 W/m- as a function of 6, one observes 
that the magnetic field at Paranal is just strong enough to col- 
lapse the populations for 6 — > n/2. In other words, by scanning 
6 across the sky with a narrow-band 20-W-class laser, we can 
observe sodium excitation all the way from full optical pumping 
to no pumping at all, leading to a severe return flux penalty, as 
shown in the following figures. 

Table 3 lists the sensitivity of tp w.rt. a 1 % perturbation in 
some simulation parameter p, more precisely .01 p) / tlr{p) - 
1, where all other parameters are those of Table 1. A value of 
b% in the last table column thus indicates that t//{p) oc p'' in 
some range around the chosen value of p. We have selected only 
those parameters that we will not study in greater detail in the 
following subsections (except /). In addition, we have excluded 
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those parameters around which ij/ is stationary (iff'ip) = 0), and 
also those parameters whose influence on the observed return 
flux is obviously linear, such as Cns- 

From Table 3, we notice that with the parameters of Table 1 
the sensitivity to changes in / is small. A comparison with 
Fig. 3b (dashed blue curve) shows that iff(I) is very shallow near 
/ = 46W/m^ (the slope is actually weakly positive due to the 
magnetic resonance bump), meaning that *?(/) - oc I, i.e. 
the absolute photon return grows linearly with the irradiance. 
The dependence of i/r on B, however, is strong: Switching from 
B = 0.23 (Paranal) to B = 0.48 (SOR, factor 2.11) decreases 
i/f by a factor of 0.76! The influence of temperature is signifi- 
cant as well since the width of the velocity distribution scales 
like T-^^ , and conversely its normalization (the number of Na 

atoms per velocity class) scales like T-^^'. However, increasing 
TNa from 185 K to 200 K, which may be a typical seasonal vari- 
ation, reduces i/r by only 3.7 %, and hence common temperature 
variations do not directly influence if/ much (we note that temper- 
ature influences the sodium abundance). Both yex and ys have a 
small influence on ^, mainly due to the presence of repumping, 
as Milonni ( I Q99) has also noted (for q - the sensitivities are 
0.049 % and 0.56 %, respectively). 

We can draw two lessons here: There is little point in 
spending much eff'ort improving the spin-exchange cross section 
'^Na-o, estimate, and, secondly, repumping makes ijj more robust 
towards some parameter variations, leading to higher stability in 
numerical simulations and possibly also in experiment. The vari- 
ation of i/r with yvcc is much larger than with ys , and at first sight 
surprisingly, the derivative is positive. We explain the positivity 
by the fact that collisions mitigate spectral hole burning due to 
recoil (neglecting recoil, the derivative is -0.089 %). Finally, we 
have included the average recoil frequency shift in the list, al- 
though it depends only on fundamental constants and A (more 
precisely hjim^^A^) - 50.0 kHz, with the atomic mass of Na 
= 3.819 X 10"^^ kg), in order to demonstrate the importance 
of proper recoil modeling at / > 1 W/m^. 

In the following subsection, we will show that the laser pa- 
rameters of Table 1 are close to optimal, given the other condi- 
tions of the table. We focus on the case - njl where the laser 
beam is directed orthogonal to the geomagnetic field (the worst 
case, for which LGS lasers must be designed to achieve a given 
return flux requirement). 



4.2. Optimization ofi// 

4.2.1. Polarization ellipticity angle 

Figure 5 plots i// asa function of the laser polarization ellipticity 
angle x- To our knowledge, this is the first publication to dis- 
cuss ijrix) for LGS. The axis limits +1 of the quantity 4-x/n de- 
note circular polarization (the poles of the Poincare sphere; in 
our convention, 4x/7t - +1 pumps towards increasing m), and 
X - indicates linear polarization. The curves pertain to differ- 
ent combinations of 6 and </>, the polar angle and azimuth of the 
B vector in a system where the laser points along z (note that for 
6 = 0, the azimuth (p becomes meaningless at any light polariza- 
tion. Conversely, when using any laser that is not purely circular 
polarized, does influence the return flux, which is often forgot- 
ten. When using purely linear polarized light, = 7r/2, = is 
equivalent to setting B = within our convention, while (f> = 7r/2 
induces the strongest effect from Larmor precession.) 

The cosine-like shape of ij/ix) presents an initially gradual 
decrease of the return flux from 4x/pi = +1, meaning that high 
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Fig. 5. i/f as a function of laser polarization ellipticity angle x^ 
where ^xln = +1 (0) marks circular (linear) states. The curves 
indicate different combinations of the magnetic field polar an- 
gle Q and azimuth angle 0. Note that i/'(y) is highest for circular 
polarization in all cases. 
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Fig. 6. Contour plot of \l/(x,0) in ph/s/sr/(W/m^) for magnetic 
azimuth = ;7r/2 at / = 46 W/m^. 



polarization purity is not required when pumping the sodium 
with circularly polarized light. If one is willing to accept a 
decrease in i/r of 5 % along the solid red curve, it is enough 
to maintain A\x\l7i > 0.8, i.e., down to an ellipticity angle of 
\x\ = 36°. The polarization extinction ratio (PER) in dB is given 
by PER = -20 logjQ tan \x\ and assumes the values of and oo 
for circular and linear polarization, respectively. The condition 
of Mx\l^ ^ 0.8 implies PER < 6.4 dB. This insensitivity is good 
news for the design of launch telescope optics. 

Figure 6 shows a contour plot of ij/ix, d) in ph/s/sr/(W/m^) 
for = 7r/2 (worst case) under standard conditions. Current LGS 
systems, even SOR's FASOR laser, operate at lower 'P/2' irra- 
diances than /p/2 = 46 W/m^ and/or at linear polarization, for 
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Fig. 8. Contour plot of if/ in ph/s/sr/(W/m-) as a function of / and 



which il/{0) varies slowly. At higher / and circular polarization, 
i/f(0) is a steeper function (along the right vertical plot edge). 
For 6 > 45°, also \l/(x) varies slowly near A\x\ln - 1, indicat- 
ing again that the circular polarization purity does not have to be 
very high. 



4.2.2. Repumping frequency offset 

Figure 7 shows as a function of A/^^, the repumping frequency 
offset between Dia and Dab, relative to 1. 7178GHz, for two 
single-frequency laser lines and ^ = 0.12at/ = 1 W/m^ (dashed 
green), / - 10 W/m- (dash-dotted blue), and 1-46 W/m^ (solid 
red). The repumping from {S,F-l) — > {S,F = 2) requires exci- 
tation to the states (P, F' - I) or (P, F' - 2) since no transitions 
(S,F^l) (P,F' = 3) are allowed. 



400 



300 



E 
o 



200 




0, circular 
Tt/2, circular 
Tt/2, circular, SOR 
(|) = Tt/2, linear 



20 30 

q{%) 



40 



50 



Fig. 9. i/r as a function of repumping power fraction q for differ- 
ent field angles 0. Dash-dotted red curve; standard conditions, 
solid blue: 0-0 (circular polarization, no Larmor precession), 
dashed black: standard conditions, but B - 0.48 G (SOR), soUd 
green: linear polarization (0 - <p - 7t/2). 



Figure 8 is a contour plot of if/(I, Afai,), zooming into the re- 
gion of Afah = 1717.8 ± 20 MHz. In order to harness the full im- 
provement of repumping, one needs to adjust Afat, with a preci- 
sion of a few MHz, particularly when at intermediate irradiance 
(/ w 10 W/m^). At larger irradiance, power-broadening washes 
out the peak (a laser beam has irradiances < / < /,nax in the 
mesosphere, located along a horizontal line in the diagram). The 
tolerance in Afat is so narrow because repumping is most effi- 
cient if a given Na atom that has been downpumped to the F=l 
lower ground state can be reexcited before its velocity changes 
due to a collision, for which the frequency offset of the repump- 
ing laser line should lie between the (S,F =1) — > (P,F' = l,2) 
transition frequencies minus the (S,F -2) {P, F' - 3) transi- 
tion frequency that the main laser line is tuned to, i.e., close to 
1717.8 MHz. 



100 4.2.3. Repumping power fraction 



Figure 9 portrays the calculated i//{q) for different magnetic field 
polar angles and field strengths B. The dash-dotted red curve 
is for standard conditions of Table 1, the solid blue curve for 
= (i.e., circular polarization and no Larmor precession), the 
dashed black curve is for standard conditions, but with the mag- 
netic field at SOR, and the solid green curve is for Unear polar- 
ization (worst case 0=<p = n/2). 

All curves peak around q = 8-25 %, and the maximum oc- 
curs at larger q the stronger the magnetic field effect is. Note 
that tfr{q) peaks at ^ > even in the absence of Larmor preces- 
sion (blue curve) since there is always a chance that atoms decay 
to the lower ground state. Moreover, iff{q) is very steep at ^ = 0, 
rewarding even the weakest repumping, which is the reason why 
broadline lasers (a; 2 GHz) often do reasonably well, despite their 
poor spectral overlap with cr(v) (see Fig. 1). Our simulations, 
however, indicate that a narrow-band laser with 12 % repumping 
beats the return flux of a single-line 2 GHz-linewidth laser by a 
factor of 3.7 under standard conditions (see next subsection). 

For circular polarization we find t//{q = 0. 12) /tlr{q = 0) ~ 3.5 
at / = 46 W/m-, underlining again that repumping is extremely 
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Fig. 10. a,b) Contour plots of q) in ph/s/sr/(W/m^) at Paranal for linear and circular polarization, respectively; d) The ratio 
i/'(circular)/i/'(linear) at Paranal; e) Same as b), but for SOR geomagnetic field; c,f) Repumping gain , q) / , q - Q) at Paranal 
and SOR, respectively, circular polarization. Blue crosses; / = 46 W/m-, q = 0.12. 



beneficial, in particular in future LGS systems involving high 
irradiance. 



Figure 10 shows six semilogarithmic contour plots of q) 
(note the different color scales of the plots). The blue crosses 
indicate the point (/ = 46W/m^, q = 0.12). Plot a) exhibits if/ 
for linear polarization at Paranal for 6 - (p - n/2, which only 
varies weakly with B. At the blue cross, if/ = 169ph/s/sr/(W/m-) 
(for q' = 0, it is only 44ph/s/sr/(W/m^)!). Plot b) shows the same 
for circular polarization {ip = 258ph/s/sr/(W/m-) at the blue 
cross, which is close to the maximum of t// in Plot b) and lies 
within the magnetic resonance 'island'). Plot d) shows the ratio 
q, 4-x/n- 1 )/(/'(/, q,x - 0). Circular polarization at 6 - n/2 
only improves the return flux at Paranal if / > lOW/m- (at 
smaller 0, there is a much higher gain, as indicated by Fig. 3). 
The improvement at the blue cross is 1 .52 (with the B-field at 
SOR, it would be only 1.2, but can reach 2.0 for smaller 6). We 
have limited the plot range in the lower right corner of d) to 1.8 
for better rendering; in reality, the value exceeds 3. Plot e) shows 
the same as b), but for the B-field of SOR. As noted in Fig. 9, the 
optimal q lies closer to 19 % under this condition. Note that even 
at q' = 0.5 (equal power in D2a and D2b), iff ^ 164ph/s/sr/(W/m-) 
is not much below its optimum of 201. Plot c) depicts the gain 
in i// when turning on repumping at standard conditions, i.e., 
if/{I, q)/t//{I, q - 0). Plot d) shows the same for SOR. At the blue 
crosses if/{q = Q.l2)/if/{q - 0) equals 3.6 and 4, respectively. 



4.2.4. Laser linewidth 

So far, we have only simulated single-frequency line lasers. We 
now turn to lasers with Lorentzian lines of finite bandwidth. 
Figure 1 1 shows two contour plots of if/{I, A/) for conditions 
at Paranal and SOR. The specific return is highest near the 
dashed blue lines, hence near 10W/m-/MHz, or 100 W/m^ per 
natural linewidth of sodium of 9.8 MHz for Paranal and near 
17 W/m^/MHz for SOR. For irradiances below about 20W/m^ 
this linearity does not hold anymore and A/ — > will yield 
higher return flux. Conversely, at / » 1 000 W/m^, A/ is no longer 
negligible compared to the Doppler width, and hence the ab- 
sorption cross section falls off in the wings of the laser spectrum 
(our simulations for q > lose accuracy for wide bandwidths, 
see the discussion of Fig. 1). In the upper left region of the plot, 
i// decreases due to low pumping that is too weak to overcome 
the population equalization from Larmor precession, and in the 
lower right we suffer from saturation. We note that d'Orgeville 
(2000) has conducted a related study (d'Orgeville's Fig. 8 differs 
from ours in that its horizontal axis shows laser power and not 
irradiance, and it is rendered in double logarithmic scale). 

Based on the dashed blue lines in Fig. 11, the opti- 
mum FWHM laser bandwidth can be estimated by A/opt ~ 
(0.025 Ip/2 - 0.5)/ B, where A/opt is in MHz, Ip/2 in W/m^ with 
20 < Ip/2 < 1000 W/m^, and B in Gauss. We expect this op- 
timal spectral power density to only depend on the relaxation 
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Fig. 11. a) Contour plot of iA(/,A/) in ph/s/srAW/m^) for 
standard conditions (Paranal). b) Same, but the geomagnetic 
field of SOR. The optimal spectral irradiance is near 10 and 
17 W/m^/MHz, respectively (dashed blue lines). 



time scales. This scaling relationship is good news for future 
LGS systems, since it means that we can achieve high return 
flux even at very large irradiances, as long as we keep increasing 
the laser bandwidth. Lasers with microsecond pulses designed 
for mesospheric spot tracking or lidar with duty cycles of 1 :20 - 
1:100, and hence high peak powers, may take advantage of this 
possibility. 
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Fig. 12. Beam efficiency Sce as a function of FWHM laser band- 
width (A/) for a P = 20 W laser at Paranal (solid red curve, stan- 
dard conditions), Mauna Kea (dashed green), and SOR (dash- 
dotted blue), depending on the local geomagnetic-field strengths 
{P - 20W, 6 - 7t/2). Dotted magenta (solid black): Laser at 
Paranal forP = 217W(P=435W). 



4.3. Optimization of laser beams 

In this subsection, we extend the analysis to entire laser beams 
by carrying out weighted integrals of i/r(/) for (Gaussian) spot 
profiles, from which we compute Sce (Eq. (7)). 

To begin, we conduct a comparison with the beam return 
flux computed by Milonni (1999). We apply our Bloch code to 
Milonni's Na conditions of B ^ 0.5G, q ^ 0, ^ 30°, ^ = 0°, 
f = 1 W, 2 X Ip/2 - 0.52 W/m^, circular polarization (4\x\/7t - 
1), no recoil, yex = 0, = 200 K, T,, = 1.0, yvcc - 0, ys - 
1/(100 /is)) and obtain Sce = 335ph/s/W/(atoms/m^), in contrast 
to Milonni's value Sce = 240ph/s/W/(atoms/m^) from Eq. (30) of 
that paper If we switch to ys = 1 / (640 /us) (long-dashed curve 
in Milonni's Fig. 6a), we obtain Sce - 359ph/s/W/(atoms/m^). 
For the combination yvcc = 1/(100/ys), ys = 1/(640/ys), we 
find See - 303 ph/s/W/(atoms/m^), underlining the importance 
of modeling diffusion in atom velocity space. We suspect that the 
discrepancy in the flux results is partly due to Milonni's (1999) 
greatly simplified S-damping formula Eq. (18), which tends to 
overestimate spin relaxation. 

Figure 12 shows a plot of the beam efficiency Sce for different 
geomagnetic field strengths at ^ = 30° (solid red: Paranal, dashed 
green: Mauna Kea, dash-dotted blue: SOR). The dotted magenta 
and solid black lines show Sce for Paranal at P = 217 W (Ip/2 = 
500 W/m^) and P = 435 W {Ip/2 = 1000 W/m^), respectively. At 
these irradiances, the ScdAf) peaks at 38 MHz and 75 MHz, re- 
spectively, without losing much of its magnitude compared to 
P = 20 W. We expect that pulsed lasers with pulse durations 
much longer than 16 ns can enjoy this high specific return. Note 
that the numerical similarity of s^e < 245ph/s/W/(atoms/m^) 
with typical values of tf) in ph/s/sr/(W/m^) is a mere coincidence, 
and the two quantities should not be confused. 

We add the caveat that the laser bandwidth is defined here 
as the short-term bandwidth, as opposed to the apparent band- 
width obtained from a long-term spectral measurement of a laser 
whose central frequency is slowly drifting, or being chirped. 



14 



R. Holzlohner et al.: Optimization of cw sodium laser guide star efficiency 



5. Conclusions 

In this article, we describe a modeling method to accurately de- 
termine the photon return flux from sodium LGS using Bloch 
(density matrix) equations. Numerical simulations are crucial, 
since it is not easy to achieve the desired optimization experi- 
mentally. 

We summarize some optimization conclusions for cw 
sodium LGS lasers as follows: 

1. The next generation of 20W-class (launched) cw laser 
sodium LGS will achieve unprecedented mesospheric irra- 
diances near /p/2 = 50W/m^ {Ipii - '50% power in the 
bucket irradiance') in good seeing conditions, using launch 
telescopes with 30-50 cm clear aperture diameter 

2. Atomic Larmor precession due to the geomagnetic field B 
can completely suppress optical pumping and thus strongly 
reduce the LGS return flux, depending on the orientation be- 
tween the laser beam and the field lines. LGS systems must 
be dimensioned for the beam pointing that yields the lowest 
return, i.e., laser beam orthogonal to B. 

3. In the future LGS irradiance regime, circular polarization 
will significantly elevate the return flux, with a rather relaxed 
polarization extinction ratio of PER < 6.4 dB. The polariza- 
tion chirality (sense) does not matter 

4. Repumping (excitation of the D2b line) with 10-20% of 
the laser power in combination with circular polarization in- 
creases the return flux further by a factor of 3^. Repumping 
can moreover alleviate spatial power-broadening and hence 
reduce the observed LGS spot size. 

5. The optimum FWHM laser bandwidth can be estimated for a 
Lorentzian line shape by A/opt ~ (0.025 /p/2 - 0.5) /B, where 
A/opt is in MHz, /p/2 in W/m^ (20 < /p/2 < 1000 W/m^), and 
B in Gauss. As a rule of thumb, the laser bandwidth in MHz 
should approximately equal the launched laser power in Watt 
divided by six, assuming a diffraction-limited spot size in 
good seeing. Employing lasers with a single line of 2-3 GHz 
bandwidth is usually several times less efficient. 

6. A laser with the above properties can achieve a specific re- 
turn of ice = 200 - 250photons/s/W/(atoms/m^) at a low- 
B location such as in northern Chile and ice - 150 - 
200photons/s/W/(atoms/m^) at a high-B location such as 
New Mexico (all values for laser beam orthogonal to B). 

We expect a numerical accuracy of +10% in our results of i// 
for a given parameter set. However, many physical simulation 
parameters such as spin relaxation and velocity-changing colli- 
sion cross sections are still to be measured more accurately (at 
mesospheric temperatures). The corresponding relaxation rates 
depend on gas pressure and temperature and thus vary exponen- 
tially with altitude. In practice, the absolute return flux uncer- 
tainty is completely dominated by sodium abundance and alti- 
tude distribution fluctuations (for instance, raising the altitude by 
6 km causes a reduction of yvcc from 1/(35 /us) to l/(100yL(s), di- 
minishing ice from 245 to 205 photons/s/W/(atoms/m^)). The ac- 
curate return flux computation of a realistic LGS system requires 
averaging along the beam, accounting for altitude-dependent 
temperature, gas concentrations, and laser irradiance. Despite 
these caveats, we do not expect the qualitative conclusions of 
this work to change when the parameter values are made more 
precise. 

We also compare the Bloch simulation results with our 
Monte Carlo rate-equation simulation Exciter and find good 
agreement. However, such Monte Carlo rate-equation methods 
are based on numerous assumptions, and they always need to be 



validated against more rigorous methods such as a Bloch equa- 
tion calculation. Moreover, our Bloch code is orders of mag- 
nitude faster than Exciter, with the runtime advantage quickly 
growing with the desired level of accuracy. Our Bloch simula- 
tion code can simulate any alkali and is publicly available on our 
websites (ESO web and Rochester online, requires Mathematica 
V. 6 or higher). 

In this work, we treat only steady-state sodium excitation 
with cw lasers. Pulsed lasers, with pulse durations close to or 
shorter than the sodium lifetime of 16 ns, will be dealt with in a 
forthcoming publication. However, we surmise that our present 
results are valid for microsecond laser pulses that are useful for 
mesospheric spot tracking and lidar On-sky experiments to val- 
idate our results with a mobile lOW-class laser unit are being 
prepared at ESO at the moment, and are expected to yield re- 
sults in 2010. 

In conclusion, our simulations indicate that we have only be- 
gun to realize the full capability of sodium LGS, to be harvested 
by upcoming generations of laser-assisted AO, both with cw and 
pulsed lasers. The future of sodium LGS, so to speak, looks very 
bright. 
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